arXiv:1508.00347vl [math.NA] 3 Aug 2015 


Manuscript 


Orthotropic rotation-free thin shell elements 

Gautam Munglani Roman Vetter Falk K. Wittel Hans J. Herrmann 


Abstract A method to simulate orthotropic behaviour in 
thin shell finite elements is proposed. The approach is based 
on the transformation of shape function derivatives, result¬ 
ing in a new orthogonal basis aligned to a specified preferred 
direction for all elements. This transformation is carried out 
solely in the undeformed state leaving minimal additional 
impact on the computational effort expended to simulate 
orthotropic materials compared to isotropic, resulting in a 
straightforward and highly efficient implementation. This 
method is implemented for rotation-free triangular shells us¬ 
ing the finite element framework built on the Kirchhoff- 
Love theory employing subdivision surfaces. The accuracy 
of this approach is demonstrated using the deformation of 
a pinched hemispherical shell (with a 18° hole) standard 
benchmark. To showcase the efficiency of this implemen¬ 
tation, the wrinkling of orthotropic sheets under shear dis¬ 
placement is analyzed. It is found that orthotropic subdi¬ 
vision shells are able to capture the wrinkling behavior of 
sheets accurately for coarse meshes without the use of an 
additional wrinkling model. 

Keywords Finite elements • Rotation-free shells • Or¬ 
thotropic materials • Subdivision surfaces • Wrinkling 


1 Introduction 


Shells possess a unique curved shape which allows them 
to carry transversal loads primarily by in-plane forces [28]. 
This feature permits them to be extremely slender if neces¬ 
sary, resulting in a wide range of potential applications. Ef¬ 
forts to understand their mechanical capabilities have there- 
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fore been crucial in the development of load carrying struc¬ 
tures. In conjunction with the finite element method (FEM), 
shells can be numerically simulated to investigate how mem- 
brane-like structures manage different loading conditions, 
thereby greatly assisting in their design {]]]. 

Biological structures like plant cell walls and skin often 
pos sess very thin sheets which undergo large deformations 
[ 10] that are difficult to capture using traditional simulation 
methods. These sheets are generally heterogenous in nature 
with anisotropic material properties |12]. Some even have 
preferential orientations due to the presence of crosslinked 
polymers, giving them orthotropic or transversely isotropic 
material properties [|3 ■ Due to their geometry, these sheets 
are often modelled with computationally efficient thin shell 
finite elements 03 . It is therefore essential that the type 
of shell elements used to approximate these structures are 
capable of simulating non-isotropic material properties. 

A traditional way to treat thin shells numerically is with 
the use of the Kirchhoff-Love theory. This theory uses a set 
of assumptions to effectively capture the properties of ge¬ 
ometrically exact thin shells in terms of its middle surface. 
These assumptions are that the shell is sufficiently thin, and 
that straight material lines initially normal to the middle sur¬ 
face of the shell remain unstretched, straight and normal af¬ 
ter deformation such that the normal stress in the thickness 
direction is negligible. Shells simulated with this approach 
have been shown to be quite accurate on a variety of numer¬ 
ical benchmarks sn. 

The choice of shape functions is crucial in the analy¬ 
sis of thin shells of Kirchhoff-Love type due to the well- 
researched C 1 continuity requirement 114011 . This requirement 
is due to the presence of second-order derivatives of dis¬ 
placement, which leads to a fourth-order equilibrium equa¬ 
tion. This in turn calls for continuous first order derivatives 
across element boundaries. In order to tackle the continuity 
requirement, higher-order conforming shape functions with 
additional degrees of freedom in Hsieh-Clough-Tocher tri- 
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angles [|7l] and Hermite quadrilaterals [40] have previously 
been used in shell elements. Alternatively, this requirement 
has been ignored and non-conforming elements with C° con¬ 
tinuity have been combined with other assumptions [29]. 
Recently, isogeometric analysis using NURBS [16], as well 
as subdivision surfaces Ig, |5[] have also been developed to 
satisfy these requirements. Subdivision surfaces in particu¬ 
lar have been implemented in triangular shells with a rotation- 
free formulation, i.e. requiring only displacement degrees of 
freedom at the mesh nodes. This rotation-free formulation 
greatly saves on computational cost due to the significant re¬ 
duction in the number of degrees of freedom, thereby gain¬ 
ing prominence in recent years lEhl 2ot 2ii, [2]. Due to their 


formulation, triangular rotation-free elements like subdivi¬ 
sion shells pose a unique set of challenges to the simulation 
of orthotropic material properties required to model deform¬ 
ing sheets B34 , 35]. 

In this article, we present a straightforward method to 
simulate orthotropy in Kirchhoff-Love rotation-free shell 
elements. This method is based on the transformation of 
the directional derivatives of the element shape functions. 
It is computationally efficient, independent of the choice of 
shape functions, and can be coupled to a wide range of prob¬ 
lems without further alterations. We use Loop subdivision 
shell elements |22] to demonstrate its accuracy for geomet¬ 
rically nonlinear problems, which include the pinched hemi¬ 
sphere benchmark and the problem of thin sheets wrinkling 
under applied shear. 

This article is divided into the following sections: first, 
the typical continuum formulation for geometrically exact 
Kirchhoff-Love thin shells is described. Section 3 then pro¬ 
vides a short overview on the development of the basis ori¬ 
entation followed by a detailed explanation of the transfor¬ 
mation method. This is followed by a few words on the con¬ 
stitutive law and discretization. Finally, the accuracy of the 
method using the pinched hemisphere benchmark is demon¬ 
strated which culminates with the application of orthotropic 
subdivision shells to the analysis of thin sheets subjected to 
shear. 


2 Kinematics 


To begin, the standard formulation of geometrically exact 
Kirchhoff-Love thin shells is revisited. The deformation of 
the thin shell is constructed using the classical stress-resultant 
formulation in curvilinear coordinates 13 ill . This formula¬ 
tion describes the shell in terms of its middle surface by in¬ 
tegrating the stress through its thickness. 

We begin by introducing some standard notation. The 
reference geometry of the shell is characterized by a mid¬ 
dle surface Q. C E 3 , with boundary F = d£l and thickness 
h which is considered to be small compared to the planar 
dimensions of the shell. With the action of applied loads, 


the shell deforms into a new configuration characterized by 
an altered middle surface £2 C E 3 . The position vectors of 
an arbitrary material point in the shell, denoted by r and r 
in the reference and deformed configurations respectively, 
are parametrized in terms of their curvilinear coordinates 
{0 1 ,0 2 , 0 3 } as 

r(0 1 ,0 2 ,0 3 )=x(0 1 ,0 2 ) + 0 3 a 3 (0 1 ,0 2 ), 
r(0 1 ,0 2 ,0 3 ) =x(0‘,0 2 ) + 0 3 a 3 (0 1 ,0 2 ), ' 

where 0 3 g [—h/2, h/2] represents the thickness coordinate, 
x and x are the parametric representations of the middle sur¬ 
faces of the configurations £2 and £2 respectively, while a 3 
and a 3 are the shell director vectors in the reference and de¬ 
formed configurations respectively. The change in the mate¬ 
rial point from x to x is described in terms of a displacement 
field u = x — x. This parametrization allows the description 
of the middle surface coordinate basis which span the tan¬ 
gent space T as 

_ 0x _ dx 

a “ = ~dQa = X ’“ ’ a “ = ~dQa ~ x -“ ’ ( 2 ) 

where Greek indices indicate integers 1 and 2, and a comma 
indicates partial differentiation. In addition, it should be noted 
that henceforth, Einstein summation is used and Latin in¬ 
dices indicate integers from 1 to 3. By applying the Kirch- 
hoff constraints, the shell directors can be described as 


ai xa 3 - _ 

a 3 = —=—, J = ||ai x a 2 1| 

a l x a 2 . I, || 

a 3 =—-—, / = ||ai x a 2 1| 


(3) 


where J is the Jacobian determinant. The infinitesimal area 
element is now expressed as df2 = 4 d 0 1 d 0 2 . By differenti¬ 
ating Eq. ©■ the covariant basis of an arbitrary point in the 
shell can be described as 


dr 

00 “ 

dr 

00 “ 


dx 

00 “ 

dx 

00 “ 



a 3 

dr 

dx 

+ u 

00“ ’ 

00 3 ' 


+ 0 3 

a 3 

dr 

dx 

00“ ’ 

00 3 

00 3 ' 


Substituting Eq.Q in the above yields 


(4) 


g« = a a + 0 3 a 3 . a , g 3 = a 3 , 

g« = a a + 0 3 a 3 . a , g 3 = a 3 , 


(5) 


while the corresponding covariant coefficients of the metric 
tensors are 


8ij g/ * g j 5 8ij g( ' g j • (ft) 

The covariant coefficients of the surface metric tensor are 
a ap = a a • , a a p = a a ■ a^ , (7) 
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while the covariant coefficients of the shape tensor are 


bal 3 = — 33 , a • &p = 83 • a a p , 

bafi = a;’>. a • ap = 83 • a a p . 


( 8 ) 


This formulation uses the Green-Lagrange strain (E/j), whose 
coefficients in curvilinear coordinates are defined as 


Eij = - {gij-gij). 


(9) 


Using additive decomposition, the Green-Lagrange strain 
tensor can be split into the non-zero coefficients of the mem¬ 
brane a a p and bending ji n p strain tensors. Neglecting higher 
order terms (<^(0 3 ) 2 ) and considering Eq.© and ©, we ar¬ 
rive at 



Fig. 1 Thin shell with reference thickness h represented by its middle 
surface £2 with non-orthogonal basis a a spanning tangent space T q at 
point q 


Eap — ®ap "h ® Pap 1 (10) 

where the coefficients of the strain tensors are 

^ap = ~^{ a ap — a ap) 1 Pap = b a p — b a p ■ (H) 

This neglection of the higher order terms makes this a first 
order shell theory, which is valid for small thickness h. 


3 Shape function derivative transformation 

Now that the basic notation has been reviewed, we describe 
how the element basis is traditionally oriented based on the 
shape functions, followed by our approach to simulate or- 
thotropy. As stated in the previous section, the tangent space 
T of the middle surface of a shell is spanned by the in-plane 
basis vectors of the reference and deformed configurations 
a« and a« respectively. When the continuum is discretized 
into elements, each quadrature point q of an element ac¬ 
quires its own tangent space which we shall refer to as T q 
as can be seen in Fig. Q. It should be noted that the number 
of quadrature points per element has no bearing on our pro¬ 
posed method as the transformation occurs on all available 
T q . Due to this property, the reduced integration approach of 
a single quadrature point per element is used henceforth. 

By construction, triangular elements generally have non- 
orthogonal basis vectors in T q . This also applies to curved 
regular meshes, where the orientation of each triangular el¬ 
ement varies from its neighbours, resulting in basis vectors 
that are not aligned between elements (Fig.©)). Isotropic 
materials do not depend on the orientation or orthonormality 
of the basis as their properties remain constant in all direc¬ 
tions, unlike the case for orthotropic materials. Therefore, to 
confer orthotropic properties to a surface, the element basis 
vectors are transformed to an orthogonal basis and aligned 
to a preferred direction across the surface as seen in Fig.©). 


The discretized form of the basis vectors can be written as a 
linear combination of shape function derivatives 


n s f n s f 

at = £ N, A x 7 , a 2 = £ N, n x/, 

7=1 7=1 


( 12 ) 


where {NiJ = l,...,n s f} are the shape functions and x/ are 
the nodal positions, with n s f being the number of mesh nodes 
accounted for by each quadrature point. These derivatives 
are defined by directions £, rj £ E 3 along the sides of the 
standard triangular master element as shown in Fig.©. Ba¬ 
sis vectors are in fact manifold versions of directional deriva¬ 
tives taken along the sides of this master element. Techni¬ 
cally, the basis can be transformed directly using an orthogo- 
nalization method like the Gram-Schmidt process 14ll . This 
would result in an orthonormal basis that spans the same 
space T q as the original shown in Fig. CQ). However, this ba¬ 
sis would be mismatched with the shape function derivatives 
shown in Eq. (fl2l) and its second derivatives. Alternatively, 
one could directly transform the covariant strain tensors a a p 
and P a p in Eg. (1111 and the contravariant resultant stresses 
shown later in Sec.© using the principal directions of or- 
thotropy @]. This procedure would need to be performed at 
every time step, which can be very computationally costly. 
Furthermore, the shape function derivatives would now be 
mismatched with the strains and resultant stresses. 

A straightforward procedure to create a basis matching 
with the first and second derivatives of its shape functions 
is to use a transformation that directly alters the first and 
second derivatives of the shape functions Ni in £2 leading 
to a new orthogonal basis without further adjustments. By 
ensuring that the new basis vectors in each T q are individ¬ 
ually aligned to the preferred direction, this method is able 
to simultaneously transform and align the basis of triangu¬ 
lar elements. An advantage of this procedure is that it only 
needs to be performed once for every element in the mesh 
for the reference configuration in a Total Lagrangian formu¬ 
lation. Other works |34j,|35] have used alternate transforma- 











4 


Gautam Munglani et al. 




Fig. 2 (a) The orientation of the basis is originally non-orthogonal and 
non-aligned, (b) The transformed basis is orthonormal and aligned to 
the preferred direction 


tion matrices built using the deformed configuration. This 
implies that, as with the direct transformation of the strains 
and resultant stresses, the process of orthogonalization is re¬ 
peated for every deformed configuration. 

The first step in this procedure is to introduce a vector to 
indicate the preferred direction of the material. This shall be 
denoted by the unit vector d € E 3 , which is not parallel to a 3 
anywhere on the middle surface £2. In order to ensure that 
the reference basis of each element is aligned, d is identical 
for every quadrature point in the mesh. The curvature of the 
shell surface predicates that d will most likely not lie in T q . 
Therefore, it needs to be projected onto T q using the relation 

&i =d-(d-a 3 )a 3 , (13) 

where ai is the new unsealed in-plane basis vector that points 
towards the preferred direction while lying in T q . The other 
new unsealed basis vector a 2 points in the direction perpen¬ 
dicular to the preferred direction, henceforth referred to as 
the perpendicular direction 

a 2 =a 3 xai. (14) 


■n V 



Fig. 3 The basis transformation corresponds to a change in the shape 
of the master triangle from £ and tj to the directions stipulated by £' 
and 7}' 


The unsealed basis a a now needs to be rescaled by the 
lengths of the original basis vectors a« to the new in-plane 
basis in the reference configuration, a' a . This is required as 
the basis transformation is non-orthogonal, thereby not pre¬ 
serving vector lengths. The angle between the original basis 
ai and the new preferentially aligned d] is given by 


0 = arctan 


ll&i x ai|| \ 

at-ai ) 


(15) 


Since J has to remain constant, a 3 can be treated as the ref¬ 
erence around which the original basis is rotated. Using the 
relation ||a 3 || = ||ai|| ||a 2 || sin 0, the new orthogonal basis a ; a 
can be defined as 

where “^- <16) 

Since li' a has been constructed, and a^ = a 3 , we can look 
at altering the shape function derivatives, N I <t and iV/.jj re¬ 
spectively. In their original state, | and 7] correspond to the 
three-dimensional Cartesian unit vectors ei and e 2 respec¬ 
tively. 

4 = [4i & 4a] = [1 0 0], 

n = lvi 02 03] = [o io]. 


These directions need to be modified to 4 ' and Tj' respec¬ 
tively. These new directions are found by constructing a trans¬ 
formation matrix T between aj and a k 


T = (a* <g> e k ) 1 (a* ®e k ) = [7jy] 


[T a p] O' 

0 1 ’ 


( 18 ) 


where e k denotes the in-plane Cartesian unit vectors and 
[Tap] I s a non-orthogonal 2x2 matrix representing an in¬ 
plane basis transformation with coefficients 



(19) 


which describes the shape of the new non-right angled mas¬ 
ter triangle as can be seen in Fig.(0. Application of [T a p] r 
modifies Nj ^ and N/ -q to derivatives with directions ^ and 
X]' respectively. This transformation also applies to the sec¬ 
ond derivatives of N/ if higher-order shape functions are 
used. This is shown by 

v , (iV/) = [r a/3 ] T v(A / ), 

H'(N,) = [T a p] T H(Nj)[T a p], 

where V'(A/) and H'(Ni) represent the modified gradient 
and Hessian matrix of the shape functions respectively. In 
tensor form, this can be rewritten as 


Ni,r\' 


1 

cs 

NJ ( 


1 

E- 

N I£'_ 


1 

(N 

njp 

_l 


N,£ 


(21) 
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The first and second derivatives of the shape functions 
have now been transformed, allowing the new orthogonal 
basis to be constructed by replacing Nj ^ and /V/.j) with N f ^ 
and N[jfi in Eq. fl 1 21 ). The altered second derivatives of the 
shape functions given in Eo. (l22l i are used to build the first 
derivatives of the basis a a p. The coefficients of the surface 
metric tensor and shape tensor defined in Eq.(|7][8]i respec¬ 
tively are then determined using this altered basis. This fi¬ 
nally results in updated membrane a a p and bending [i n p 
strain tensors shown in Eq.tQj}. Once these geometric quan¬ 
tities have been built, the constitutive law in curvilinear co¬ 
ordinates needs to be defined for the orthotropic case. This 
combined with the assembly of the stress resultants with the 
new orthogonal and aligned basis is described in the next 
section. 


4 Constitutive model and discretization 


The material properties of the shell are determined by the St. 
Venant-Kirchhoff constitutive law. The strain energy den¬ 
sity function for the orthotropic case is altered from the Koi- 
ter energy density functional [ 19] to 


W = ■ 


K a P 


hH a Wa a pCC r s + ^H a Wp ali li rS 


(23) 


where the coefficients of the elasticity tensor H a ^ s are 
H*W = Vi + (1 ~ Vl) (a ay a^ + d aS d^) , (24) 

and the coefficients of the stiffness K a ?‘ are 


[/if 11 K 22 K n ] 


E\ E 2 G\2 

1 — Vi V2 1 — Vi V2 1 — Vi 

(25) 


The above relation includes the elastic modulus E, the Pois¬ 
son ratio v and the shear modulus G, where the index 1 rep¬ 
resents the preferred direction d of the material, and 2 rep¬ 
resents the perpendicular direction. The moduli and Pois¬ 
son ratios are related by E\ V 2 = £2 V[ ■ The derivative of the 
strain energy density functional with respect to the mem¬ 
brane strains a a p and bending strains j 8 a p gives the resul¬ 
tant membrane stresses n a P and bending stresses m a P of the 


element respectively 
n al} = = hK a V H a ^ s a vS , 

y 3 (26) 

m aP = dW = h K<x pHaPyS^ s 

dpap 12 

Now that the form of the stress resultants has been eluci¬ 

dated, we can proceed to the discretized form of the equilib¬ 
rium equations. This is achieved by approximating the min¬ 
imization of the total potential energy, which is obtained by 
summing the contribution of the internal elastic energy (0 lnt ) 
with the external energy (0 ext ) according to 

<Mu] = <r t [u] + <T t [u], (27) 


The internal energy is the Koiter strain energy density is in¬ 
tegrated over the reference middle surface, and the external 
energy is the sum of the external load q per unit surface area 
and the traction N per unit edge length, given as 


f nt [ u]= f_Wd£2. 

Jn 

0 ext [u] = - /_q ud?T- f_ N-udf, 
Jn Jr 


(28) 


respectively. To solve the elastic energy minimization prob¬ 
lem for the displacement field u, the first variation of (j) is 
taken and augmented with the inertial term containing the 
mass matrix. This results in 


o = ff-fft + £M„ii 7) 


(29) 


where 


ff =~L y d ~y +'»“ p yy <&> 

Jn V OU; J 

f, u = l_ q/V/dl2+ /_N/V/dr. 

Jn Jr 

M u = J hpNiNjdQ. 


(30) 


In the above relation, Mjj is the mass matrix for dynamic 
analysis. Eo.(l29l> can now be evaluated element-wise using a 
quadrature rule. The constant-average acceleration method 
is used for time integration. This is obtained by setting y = 
1/2 and /? = 1/4 in the widely used Newmark family of 
methods] 23, 29j]. Further details on this setup can be found 
in Refs.|6^5[|. 

Methods that directly transform the components of f nt 
will require the additional transformation of f™ 1 to the new 
local coordinate system in order to implement boundary con¬ 
ditions which involve shape function derivatives. All rotation- 
free shells require this supplementary step for their imple¬ 


mentation of certain natural boundary conditions 1251 13 
HESS Tlil , [24] I including frictional contact [BjJ , in-plane 


shear traction based conditions [21] and displacement-dependent 







































6 


Gautam Munglani et al. 


pressure loads [30]. Some widely used formulations addi¬ 
tionally require these shape function derivatives for curva¬ 
ture gradient calculations [25, 20] or hourglass stabiliza¬ 
tion I2]]. Our proposed method does not require these further 
transformations. 

The details of the simulation method and its advantages 
have been presented and are now followed by two numerical 
studies in the next section to demonstrate the efficiency and 
accuracy of orthotropic subdivision shells. 


5 Numerical studies 

We present two examples to demonstrate the accuracy of 
the described transformation method. The first is a standard 
numerical benchmark of orthotropic behaviour for geomet¬ 
rically nonlinear problems. The second showcases the effi¬ 
ciency of orthotropic subdivison shell elements by analyz¬ 
ing wrinkling behavior of sheets. Both examples have been 
simulated using a single quadrature point per element ( re¬ 
duced integration ), located at its barycenter. It has been ar¬ 
gued in Refs. AS that a single quadrature point in adequate 
for the simulation of geometrically nonlinear large deforma¬ 
tion problems with subdivison shells. 


with large rigid body rotations. While this setup was origi¬ 
nally developed for isotropic shells, it was modified to take 
into account the orthotropic case in Refs.[34, 35]. The prob¬ 
lem consists of a hemisphere with a 18° hole at its north 
pole loaded by four equal point loads on its equator 90° from 
each other. These forces are diametrically opposite in direc¬ 
tion, with a pair of tensile and compressive loads. The shell 
radius R is 10, with a thickness h of 0.04, elastic modulus in 
the circumferential direction E c of 6.825 x 10 7 , and a Pois¬ 
son ratio v c of 0.3. The elastic modulus in the meridional 
direction E m and shear modulus G cm are given in Tab. 0 
for different material properties. For clarity, the degree of 
orthotropy is defined as A = E m /E c . The meshes were con- 


Table 1 Material properties for the pinched hemisphere benchmark 


A 

E m 

^cm 

1.0 

6.825 x 10 7 

2.625 x 10 7 

0.9 

6.143 x 10 7 

2.518 x 10 7 

0.5 

3.413 x 10 7 

1.896 x 10 7 

0.1 

6.825 x 10 6 

5.884 x 10 6 


5.1 Pinched hemispherical shell 



Fig. 4 Control mesh for deformed hemisphere (A = 0.5) with its as¬ 
sociated bending energy density after being loaded at A and B with 
pairwise opposite point loads 


structed by subdividing each triangular element into four un¬ 
til a fine enough discretization was obtained. Further data 
on the relevent isotropic case can be found in Refs. |3„2. 331. 
Due to the small number of degrees of freedom (DOF) re¬ 
quired for an accurate result and the lack of symmetry, the 
entire hemisphere was modelled. The hemispherical mesh 
consists of 16 x 64 subdivision elements as used by Ref. |3il 
with 3264 DOF (3 DOF per node). Fig.© illustrates two of 
the four diametrically opposite point loads (at A and B) on 
which an increasing load with a maximum of 100 is placed. 
Tab.© shows the maximum displacements of each material 
property scenario. The orthotropic subdivision shells (SD3R) 
are compared with quadratic shell elements with reduced 
integration (S8R) from the commercially available Abaqus 
software |@] and the rotation-free Basic Shell Triangle (BST) 

134)!, [351] in Fig-©- Very good agreement with existing data 
leads us to conclude that the orthotropic subdivision shell is 
able to accurately model problems with large deformations 
and both stretching and bending with a substantially lower 
number of DOF than the BST and S8R elements found in 
literature. 


Table 2 Displacements of A and B at maximum load 100 


A 1.0 0.9 0.5 0.1 

A 5.918 6.125 7.019 8.716 

B 3.350 3.407 3.629 3.978 


This example is used to study how the shell performs 
when subjected to coupled stretching and bending stresses 
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Fig. 5 Load-displacement curves for the deforming hemispherical shell with different degrees of orthotropy A. The hemispherical SD3R mesh 
contains 3264 DOF for the entire hemisphere. For comparison with other methods, the number of DOF for a quarter of the hemisphere is (16 x 
17) x3 = 867 where 16 is the number of elements along the meridional axis. The S8R results were obtained using Abaqus (8j] and the BST results 
from Refs.j3.lH] 


5.2 Wrinkling of orthotropic sheets 


The wrinkling of sheets has been studied extensively with 
shell and membrane elements |27, 38, 17, 11]. The mem¬ 
brane element formulation needs to be augmented with a 
wrinkling model, but can be applied on coarser meshes, mak¬ 
ing it relatively computationally inexpensive. Subdivision 
shells have been used in the past to simulate the deforma¬ 
tions of thin membranes B37II . To demonstrate the ability 
of orthotropic subdivision shell elements to reproduce the 
wrinkling phenomenon on coarse meshes, the reference shear 
test described in Ref. [27] is simulated. 


The test consists of a prestressed sheet sheared by dis¬ 
placement control. The objective is to find the critical shear 
displacement ( u c ); a bifurcation point signaling the onset of 
wrinkling, and the maximum amplitude of the wrinkles at 
the maximum specified shear displacement. The 200 mm x 
100 mm sheet is first prestressed with a displacement of 1 
mm along its short axis, and then sheared by continuously 


displacing its upper edge by 10 mm along its long axis. The 
lower edge of the sheet remains pinned for the duration of 
the simulation (Fig.©). The material is 0.2 mm thick and 
has a preferred orientation of a = 30° to the v-axis. The ma¬ 
terial properties for the isotropic case are E = 600 N/mm 2 
and v = 0.45 [27[], while for the orthotropic case E\ = Eo 


2 L,\ = 


Meshes 


106.6 N/mm 2 , v = 0.22 and G n = 11 -3 N/mm 2 | 
with resolution ranging from 8 x 4 (135 DOF) to 56 x 28 
(4959 DOF) are simulated. 


In these simulations, a wrinkle is defined as a single fold 
in the sheet, which begins from the baseline u z = 0, rises to 
its maximum value and back to u z = 0. As seen from the 
results in Fig.©, the isotropic case shows a critical shear 
displacement of u c = 1.67 mm, with a maximum wrinkle 
amplitude at a shear displacement of u y = 10 mm of 1.9 
mm which exactly corresponds to the value obtained with 
3D membrane elements in Ref. ^27h. The number of wrinkles 
(14) for the maximum resolution of 56 x 28 elements (4959 
DOF) also agrees with the existing data. In the orthotropic 
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Shear 

u y = 10 mm 



Prestress 
u x = 1mm 


Fig. 6 Boundary and loading conditions of sheet to investigate wrin¬ 
kling behavior 127 1 


case, the critical shear displacement is u c = 1.81 mm which 
compares favorably with Ref. El . The maximum amplitude 
at a shear displacement of u y = 10 mm is 2.2 mm, which is 
slightly larger than the 2 mm specified in the original study, 
while the sheet develops 10 wrinkles following the same 
profile. 




Fig. 8 Wrinkling of sheet with (a) isotropic f27ll and (b) orthotropic 
111 properties using a resolution of 56 x 28 elements with 4959 DOF. 
The displacement in the normal direction ( u z ) is given in mm and the 
wrinkles are marked with black dashes 


Strikingly, even a coarse mesh resolution of 16 x 8 ele¬ 
ments (459 DOF) yields critical shear displacements of 1.68 
mm and 1.81 mm for the isotropic and orthotropic case re¬ 
spectively, which are very small deviations from the high 
resolution meshes. The number of wrinkles converge more 
slowly for the isotropic case, requiring the 48 x 24 (3675 
DOF) resolution to reach 14 wrinkles, while 10 wrinkles are 
already obtained for the 32 x 16 (1683 DOF) resolution in 
the ortho tropic case. 
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Fig. 7 Convergence of the critical shear displacement and number of 
wrinkles with mesh resolution 



6 Conclusions 

A method to simulate orthotropy in rotation-free shell fi¬ 
nite elements was proposed. This approach transforms the 
derivatives of the shape functions to orthogonalize and align 
the basis for each triangular element. The implementation 
was performed using Kirchhoff-Love type subdivision shells. 
This approach requires negligible computational overhead 
as the transformation is only performed once in the unde¬ 
formed configuration. The standard pinched hemispherical 
shell with 18° hole benchmark for orthotropic behaviour 
proved that the method is accurate. 

The combination of orthotropic material behaviour and 
the efficiency and robustness of subdivision finite elements 
has many advantages in simulating smooth membranes with 
arbitrary topologies. This implementation was then used to 
simulate the wrinkling of sheared orthotropic sheets where 
its efficiency was compared to existing methods. 

The shape function derivative transformation applied to 
traditional subdivision shell finite elements has been shown 
to be very well suited for the simulation of material aniso¬ 
tropies. The analysis of the sheared membrane showed that 
the wrinkling of an orthotropic sheet was accurate even for 
coarse meshes, with the added advantage that no specific 
wrinkling models are required to compensate for the lack of 
bending stiffness in membrane elements. 
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